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Abstract 

The Faddeev equations for the three body bound state are solved directly 
as three dimensional integral equation without employing partial wave de- 
composition. The numerical stability of the algorithm is demonstrated. The 
three body binding energy is calculated for Malfliet-Tjon type potentials and 
compared with results obtained from calculations based on partial wave de- 
composition. The full three body wave function is calculated as function of the 
vector Jacobi momenta. It is shown that it satisfies the Schrodinger equation 
with high accuracy. The properties of the full wave function are displayed and 
compared to the ones of the corresponding wave functions obtained as finite 
sum of partial wave components. The agreement between the two approaches 
is essentially perfect in all respects. 
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I. INTRODUCTION 



Three nucleon bound state calculations are traditionally carried out by solving Faddeev 
equations in a partial wave truncated basis, which leads to a set of a finite number of coupled 
equations in two variables for the amplitude. This is performed either in momentum space 
|0-0], in configuration space ||||, or in a hybrid fashion using both spaces ||. Though a 
few partial waves often provide qualitative insight, modern three nucleon bound state calcu- 
lations need 34 or more different isospin, spin and orbital angular momentum combinations 
@. It appears therefore natural to avoid a partial wave representation completely and work 
directly with vector variables. This is common practice in bound state calculations of few 
nucleon systems based on variational Jf| and Green's function Monte Carlo (GFMC) meth- 
ods H . Those methods are normally carried out in configuration space, and while providing 
accurate results for the binding energy, the GFMC method samples the wave function only 
stochastically. Up to now the Faddeev equations have been solved applying vector variables 
only in configuration space for pure Coulomb bound state problems, the e~e~e + and ppfi~ 
systems ||. There are other accurate techniques for solving for the three body binding 
energy based on correlated Hyperspherical Harmonic variational techniques pO|JTl|] and the 
Gaussian basis coupled channel methods |12|, which, however, all work with partial wave 



decomposition. 

Our aim is to work directly with vector variables in the Faddeev scheme in momentum 
space. Under these circumstances the two-body t-matrices should also be determined di- 
rectly as function of momentum vectors. This is not too difficult and has been recently 
demonstrated for two-body t-matrices based on Malfliet-Tjon type potentials The 



choice of momentum vectors as adequate variables is also suggested by the NN force. Here 
the dependence on momentum vectors can be rather simple, e.g., in the case of the widely 
used one-boson exchange force, whereas the partial wave representation of this force leads 
to rather complicated expressions [Q. 

Instead of first summing up two-body interactions to infinite order into two-body t- 
matrices and obtain connected kernel equations, the Schrodinger equation for the three 
nucleon (3N) bound state can be solved directly ]T3| . In that work the Schrodinger equation 
was solved in momentum space for identical particles without partial wave decomposition 
for Malfliet-Tjon type interactions. We do not want to follow this avenue but rather solve 
the Faddeev equations in momentum space directly in a three dimensional form. Section II 
describes our solution of the Faddeev equations as function of the momentum vectors and 
contrast this new approach with the standard form given in a partial wave decomposition. 
We work with spinless nucleons and thus formulate the Faddeev equations for a system of 
three bosons. In Section III we discuss details of our algorithm for solving the Faddeev equa- 
tions for the three body bound state and demonstrate the stability of our numerical approach 
as well as compare to results obtained from partial wave based calculations. In Section IV 
we show properties of the Faddeev amplitude and the full wave function and demonstrate 
that the results obtained from the three dimensional approach and the ones coming from 
partial waves agree very well with each other. We compare momentum distributions and 
expectation values calculated within the two schemes. In Section V we demonstrate that the 
Schrodinger equation is fulfilled pointwise by the full wave function generated in the three 
dimensional approach. To the best of our knowledge this has never been demonstrated 
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before with this accuracy. We conclude in Section VI. 



II. THREE-BODY BOUND STATE EQUATION 

The bound state of three identical particles which interact via pairwise forces Vi = Vjk 
(ijk = 123 and cyclic permutations thereof) is given by a Schrodinger equation which read 
in integral form 

|*) = G £>|*). (2.1) 

i=i 

Here the free propagator is given by Go = {E—Hq)^ 1 , and H stands for the free Hamiltonian. 

3 

Introducing Faddeev components |\&) = with \ipi) = GoVi\^f) leads to the three 

i=i 

coupled integral equations 

m = G t t J2\^)- (2-2) 

The operator describes the two-body t-matrix in the subsystem jk. If we consider identical 
particles (here bosons, since we are omitting spin), the three nucleon wave function has to 
be totally symmetric. As a consequence, \if>i), and \ip3) are identical in their functional 
form, only the particles are permuted. Thus it is sufficient to consider only one component 

\1>) = G tP\il>), (2.3) 

where the arbitrarily chosen index 1 is dropped. In this case the permutation operator P is 
given as P = P12-P23 + P 13^23- The complete three nucleon wave function is then given as 

|*> = (1 + P)|V>. (2.4) 

In order to solve Eq. ( |2.3| ) in momentum space we introduce standard Jacobi momenta 

p i = |(k i -k fc ) (2.5) 
q; = |(ki - + k fc )), 

where ijk = 123 and cyclic permutations thereof. With the arbitrary choice % — 1 and 
omitting the index 1 for clarity of presentation Eq. (|2.3| ) reads 

<pq|^> = rp 1 2 — 5-2 (pq|^|-0> (2.6) 

with 

(pq\tP\7p) = J dY d 3 p' d\" d 3 p" (pq|t|p'q'> (pVI^IpV) (pVl^)- ( 2 -7) 

The momentum states are normalized according to (p'q'|pq) = 5 3 (p' — p) 5 3 (q' — q). The 
Jacobi momenta in the different two-body subsystems (13) and (12) are expressed through 
those defined in the subsystem (23) via 
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qi = -P2 - |qa 

pi = -5P2 + f q2 
qi = Pa - |qs 



pi = -5P3 - |qs- (2-8) 



This allows to evaluate the permutation operator given in Eq. ( |2.7| ) as 

(p q \P\v q ) = (p q Ip q )2 + (p q |p q )3 

= <* 3 (p' + |q' + q") 5 3 (p" - q' - |q") 

+ ^ 3 (P' - §q' - q") W + q' + |q")> ( 2 -9) 

where the indices 2 and 3 indicate the corresponding subsystem. Inserting this relation into 
Eq. ( |2.7| ) and then into Eq. ( |2.6| ) leads to the expression for the Faddeev amplitude which 
serves as starting point for our numerical calculations, 

<P^> = E _±l_J_2 / d V UP, |q + q', E - £q 2 ) (q + ±q', q'|V>- (2.10) 

m^ 4m * 

Here £ s (p, q, E) is the symmetrized two-nucleon t-matrix, 

* B (p, q, E) = t(p, q, E) + t(-p, q, £?). (2.11) 

We would like to mention that the so obtained Faddeev amplitude fulfills the symmetry 
relation ^(pq) = pq), as can be seen from Eq. (|2.10|) . 

In Eq. ( |2.10| ) the Faddeev amplitude ?/>(pq) is given as function of vector Jacobi momenta 
and obtained as solution of a three dimensional integral equation. For the ground state 
?/>(pq) is a scalar and thus depends only on the magnitudes of p and q and the angle 
between those vectors. In order to solve this equation directly without inserting partial 
wave projections, we have to define a coordinate system. We choose the vector q parallel 
to the z-axis and express the remaining vectors p and q' with respect to q. As variables we 
have the magnitudes of the vectors as well as the following angle relations, 

x = p • q = cos 9 
x = q' ■ q = cos 6' 

y = p-q' = cos7, (2.12) 

where 



cos 7 = cos 6 cos 6' + sin 6 sin 6' cos(<p — ip') = xx' + yl — x 2 yl — x' 2 cos ip' . (2.13) 

To arrive at the last relation we took advantage of the freedom of choice for one azimuthal 
angle and set ip = 0. We have this freedom due to the ip' integration over the full 2n interval. 
With these variables Eq. (p. 10 ) can be explicitly written as 

ip{p,q,x) = — — , 2 — _3_ J dq'q' 2 J dx' J dip' 

m" Am" q _i q 

X ts (p, vV + q' 2 + qq'x', E - ^ 

1 „> s 



X i, Lq 2 + \q' 2 + qq'x', q' , ] , (2. 14) 

V |q+ ?q 
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where ||q+ q'| = \J\q 2 + q' 2 + qq'x' and |q+ |q'| = \Jq 2 + \q' 2 + qq'x'. 

The above equation, Eq. (|2.14j ), is the starting point for our numerical algorithms, and the 

details will be described in the next Section. 

In a standard partial wave representation |i| Eq. (|2.14|) is replaced by a set of coupled 
two-dimensional integral equations 

Mp,q) 



m r Am ^ 
oo 1 



X 



Y. [*tf f<hf i ' (P ' Tl ' £ ,-^ 2) G,( q ,M %A (2.15) 

' 0-1 



= v q + iQ + qq x 



where 



-*2 = \]q 2 + \q' 2 + qq'x'. (2.16) 

Here Gu'(q,q',x f ) together with the shifted arguments tti and 7r 2 are the partial wave rep- 
resentation of the permutation operator. The explicit form of Gw(q,q' ,x') can be found in 



Ref. ||16|| . However, the expression given in this reference is more complicated due to spin 
and isospin variables. For the convenience of the reader we give Gu>(q,q',x') defined for 
general total orbital angular momentum in Appendix A explicitly. The quantities t\ are the 
partial wave projected two-body t-matrices. Due to the symmetry requirement of the Fad- 
deev amplitude, I runs over even values only. The infinite set of coupled integral equations 
given in Eq. ( |2.15|) is truncated in actual calculations at a sufficiently high value of I. The 
full Faddeev amplitude ip(p,q,x) reads then 



ip{p,q,x) = J2—^—Pi( x )Mp>Q)- ( 2 - 17 ) 



From the Faddeev amplitude q) the three nucleon wave function is obtained by 
adding the components defined in the different subsystems as given in Eq. ( |2.4| ). After 
evaluating the permutation operator P, the wave function is given as 

* (p, q) = V»(p, q) + VK-fq - §p, p - f q) + V>(|q - §p, -p - |q)- ( 2 - 18 ) 

Already here we see that \l/(p, q) = p, q), provided the Faddeev components fulfill this 
symmetry. Again the momentum q is chosen in the direction of the z-axis and after some 
algebra the explicit expression for the three nucleon wave function reads 



V(p,q,x) = ip(p,q,x) 



3„2 1^2 1 , 



+ ^ ( \ J\q 2 + P 2 + 3pqx, Jp 2 + \q 2 - pqx, 2 \ ^ 

V V V 4 | _| q _ ip||p- iq| 



+ if, (yiq 2 +p 2 -3pqx, ^ + \q 2 +pqx, ^ ± \f ) , (2.19) 

V If q — ipII - p - |q|y 
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where the magnitudes in the denominators of the angle variables are given by 

|q±fp 



f <f + P 2 ± 3pqx 
T3 



2 V 4 



p ± = Jp 2 + \q 2 ± pqx. 



The wave function is normalized according to 



which reads explicitly 



Jd 3 p d 3 q ^ 2 (p,q) = l, 



OO OO 1 

8ir 2 J p 2 dp J q 2 dq J dx ^ 2 (p,q,x) = 1. 



(2.20) 



(2.21) 



(2.22) 



The properties of the three nucleon wave function will be studied in Section IV. 
In a partial wave representation Eq. (2.18) takes the form 



*(p, q,x) = J2 1 Pi(x) *i(P> ?)» 



where 



(2.23) 



Vi(p,q) =Mp,Q) + J2 / dxG ul (p,q,x) ^(71^712) 



-1 



and 



7*1 

TT 2 



\p 2 + f^q 2 + hpx 



p 2 + \q 2 — qp x - 



(2.24) 



(2.25) 



The quantity Gw(p,q,x), which results from applying the permutation operator, can be 
found in Ref. but should be considered without spin and isospin factors. The resulting 
expression for the here considered bosonic case is given in Appendix A. 



III. CALCULATION OF THE THREE-BODY BINDING ENERGY 



For our model calculations we use Yukawa interactions of Malfhet-Tjon [18] type 



w - p) = h {w=WT& - W^WTa) ' (3 ' 1) 

We study two different types of pairwise forces, a purely attractive Yukawa interaction and 
a superposition of a short-ranged repulsive and a long-ranged attractive Yukawa interaction. 
In order to be able to compare our calculations with results obtained by other methods we 



use the Malfhet-Tjon potentials MT-IV and MT-V [[Tj|]. However, we use slightly different 
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parameters as given in Ref. [18 to compare with Refs. [§,|l!J as well as 0] and give our 
parameters in Table I. The calculated values for the deuteron binding energy are = 
-2.2087 MeV for the MT-IV potential and E d = -0.3500 MeV for the MT-V potential, 
respectively. We also need to point out that we calculate the potentials as functions of 
vector momenta and thus define the interaction as a truly local force acting in all partial 
waves. 

With this interaction we first solve the Lippmann-Schwinger equation for the fully-off- 
shell two-nucleon t-matrix directly as function of the vector variables as described in detail 
in Ref. \T3\. The so obtained t-matrix is then symmetrized to get t s (p',p,x; E — ^g 2 )- 
We would like to point out that after having solved the Lippmann-Schwinger equation on 
Gaussian grids for p, p', and x, we solve the integral equation again to obtain the t-matrix at 
points x = ±1. Thus, when iterating Eq. ( |2.14 ), we do not have to extrapolate numerically 



to angle points x of t s (p',p, x; E — ^g 2 ), which can very well be located outside the upper 
or lower boundary of the Gaussian grid of the i-matrix. 

The eigenvalue equation, Eq. (|2.14j ), is solved by iteration. Schematically this can be 
written as 

\il) = K,{E)ij), (3.2) 

and we search for E such that A = 1. The functional behavior of K{E) is determined by 
the two-body t-matrix, and A = 1 is always the largest positive eigenvalue regardless of the 
two-nucleon potential being employed. For this reason the simple iteration starting with an 
essentially arbitrary vector is sufficient. We start with a vector if)(p, q,x) ~ 1/ ((l+p 2 )(l+g 2 )) 
and stop the iteration when Eq. ( |3.2|) is fulfilled with a relative accuracy of 10~ 10 at each 
point (p, q, x). 

In order to solve the eigenvalue equation, Eq. ( j2.14|) , for the Faddeev component ip(p, q, x) 



we use Gaussian grid points in p, q, and x, though only the g-variable is a true integration 
variable. Typical grid sizes are 97 x 97 x 42 to obtain an accuracy in the binding energy of 
4 significant figures. This grid contains x points for x = ±1 as well as the points p = and 
q = 0. The iteration of Eq. ( |2.14|) requires a two-dimensional interpolation on ip(p, q, x) in the 



variables p and x. By including the additional grid points, we avoid the extrapolation outside 
the Gaussian grid. The g'-integration in Eq. ( |2.14| ) is cut off at a value of g max = 20 fm -1 . 



The integration interval is divided into two parts, (0, go)U(go, <Zmax), in which we use Gaussian 
quadrature with NQ1 and NQ2 points respectively. The value for go is chosen to be 10 fm -1 . 
Typical values for NQ1 and NQ2 are 64 and 32 to obtain the above mentioned accuracy. For 
the distribution of quadrature points we use the maps given in Ref. ptj . The x' integration 



requires typically at least 32 integration points, while for the if' integration 20 points are 
already sufficient. The if' integration acts only on t s and can thus be carried out before 
starting to solve Eq. ( |2.14| ). The p variable is also defined in an interval (0,p ) U G°o>Pmax), 



where po is chosen to be 9 fm and p max = 60 fm . The two intervals contain NP1 and 
NP2 points, and we choose NP1 = NQ1 and NP2 = NQ2 respectively. 

The momentum dependencies given in Eq. (|2.14j) suggest that we solve the two-body 
i-matrix t s (p,p', x\ e) on the momentum grid p and p' for the energies e = E — ^g 2 dictated 
by the same g-grid. It turns out that it is sufficient to choose for the variable x the same 
grid which is used for solving the two-body Lippmann-Schwinger equation for t s (p,p',x;e). 
In order to obtain the second momentum and the angle for t s (p,p',x;e) required in the 
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integration of Eq. (|2.14|) , we have to carry out two-dimensional interpolations. We use the 



Cubic Hermitean splines of Ref. pl |. The functional form of those splines is described 



in detail in Appendix B of this reference and shall not be repeated here. We find these 
splines superior to standard B-splines |22| in capturing the peak structure of the two-body 
t-matrix, which occurs for off-shell momenta p ~ p' . In Fig. 1 the off-shell structure of the 
symmetrized t-matrix t s (p,p', x; e) for a center-of-mass (cm.) energy e = —50 MeV and an 
off-shell momentum p' = 0.75 fm -1 is displayed. The two steep dips at x = ±1 dominate the 
off-shell structure and are present for all off-shell momenta p ~ p'. An additional advantage 
of the Cubic Hermitean splines is their computational speed, which is an important factor, 
since the integral in Eq. ( |2.14| ) requires a very large number of interpolations. Finally, 
to obtain the Faddeev component entering the integration we also need a two-dimensional 
interpolation and use the same as for the t-matrix. 

In Table II we show the convergence of the three nucleon binding energy as function 
of the number of grid points for the Malfliet-Tjon potentials MT-V and MT-IV. We use 
the potential parameters as given in Table I. In Table III the convergence of the energy 
eigenvalues for the same potentials resulting from the solution of the partial wave Faddeev 
equation for an optimized set of p and q grid points as function of the angular momentum I 
is given. 

As demonstrated in Table II, the calculation of the three nucleon binding energy using 
the Malfliet-Tjon potential V converges to a value of E t = —7.7365 MeV. The solution of 
the Faddeev equation in partial waves converges to E t = —7.7366 MeV. Here convergence is 
reached for / = 12. Both of the here calculated values for the potential MT-V are in excellent 
agreement with the value E t = —7.7366 MeV given in Ref. ]nj. The agreement with the 
slightly older calculation of Ref. Q , which uses 6 partial waves for the potential and yielding 
E t = —7.736 MeV is also excellent. The agreement with the value E t = —7.7365 MeV 
obtained from a the paired potential basis calculation [10] gives the same excellent agreement. 

For the purely attractive potential MT-IV our calculations of the binding energy con- 
verge to E t = —25.050 MeV as shown in Table II. The solution based on partial wave 
decomposition converges to E t = —25.057 MeV for I = 12. Both values are also in very 



good agreement with the value E t = —25.05 MeV quoted in Ref. fT5]| . As we can see from 
these comparisons to calculations of the three nucleon binding energy based on the solu- 
tion of the Faddeev equations in partial waves, our results provide the same accuracy while 
the numerical procedures are actually easier to implement. In the three dimensional case 
there is only one single equation to be solved, whereas in the partial wave case one has a 
set of coupled equations with a kernel containing relatively complicated geometrical expres- 
sions Gu'(p,p' ,x) and ti(p,p',x) matrices to larger values of /. The latter ones are driven 
by the partial wave projected potential matrix elements Vi(p, p'), which require for large 
/-values and larger p-values (80-100 fm _1 ) great care to be generated reliably. In contrast, 
the t-matrices are generated more easily as functions of vector momenta from the original 
potential. Once an accurate and fast two dimensional interpolation scheme is established, 
the implementation of Eq. (|2.14|) is rather simple. 



IV. THE THREE NUCLEON WAVE FUNCTION 



S 



A. Properties 



When solving Eq. (|2.14j) for the bound state, we obtain the Faddeev component 



ijj(p, q, x), and the symmetry relation for the Faddeev component reads explicitly ip(p, q, x) = 
ijj(p,q, —x). Note that in Eq. ( |2.14|) this symmetry is only fulfilled because of the (p- 



integration over the complete 2ir interval. We verified that our numerical solution satisfies 
this symmetry with high accuracy. In fact, this symmetry property can be implemented in 
the iteration of Eq. ( |2.14p to cut down the size of the field ijj(p, q, x) and thus save time and 



memory when computing the integral. We would like to remark that for the initial calcu- 
lations we did not take advantage of this symmetry property in order to use the numerical 
verification as a test for the accuracy of the integration and interpolation. 

In Fig. 2 we display magnitude of the Faddeev component at a fixed angle, ip(p,q,x = 1) 
obtained from a calculation based on the MT-V potential. The norm is chosen such that 
(^l 1 !/) = 3 (ip\^>) = 1, where \1/ is the full 3N wave function as given in Eq. (2.18). The 
major contributions to tp(p,q,x) arise from momenta p and q less than 1.5 fm -1 . It also 
turns out that the dependence of ip(p, q, x) on the angle between the Jacobi momenta p and 
q is so weak that it can only be detected on a logarithmic scale but not when comparing 
linear plots at different values of x. This is understandable, since the admixture of the I > 
partial waves into the Faddeev amplitude is less than 0.01%. To show this we choose the 
normalization differently from above as 



1 oo oo 

J dx J dpp 2 J dqq 2 ip 2 (p,q,x) (4.1) 

-10 



or, according to Eq. ( |2.17| ) 



oo oo 



1 = 12 J d PP 2 J dqq 2 ipf(p,q) 

1 



= J2 F i- ( 4 -2) 
i 

Similar we define for the full wave function 

oo oo 

1 = E/ d PP 2 J dqq 2 #?(p,g) 

1 

= Y,Wi- (4.3) 
i 

The partial wave contributions F\ and Wi are given in Table IV for the different I values. 

The logarithmic presentation in Fig. 2 allows to see not only the smooth bell-like decrease 
of the Faddeev amplitude for p and q less than 1.5 fm -1 but also the deep valleys caused 
by the node lines of the amplitude. A more detailed inspection of the Faddeev amplitude is 
presented in Fig. 3. Here the partial wave amplitudes ipi(p, q) are projected out numerically 
from the full amplitude ip(p,q,x) and compared to the amplitudes ipi(p,q) directly deter- 
mined from the coupled set of equations given in Eq. ( |2.15| ). We show contour plots for the I 
values 0, 2, and 4. One can see that for the higher I values the maximum of the partial wave 
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amplitudes moves to larger p and q values. It is also apparent from the figure that the graphs 
representing the three dimensional calculation and those for ones based on partial waves are 
nearly everywhere indistinguishable. This underlines the accuracy and consistency in the 
numerical realization of both schemes. Only for I = 4 and very small g-values deviations 
occur. There the contour lines have to be almost parallel to the p-axis, indicating the proper 
threshold behavior for very small q. This is not the case in the 3D approach. However, the 
value of the amplitude is smaller by 6 orders of magnitude compared to the value of the the 
I = amplitude at these values. In both cases we hit numerical inaccuracies here for the 
number of grid points chosen. 

The evaluation of the 3N wave function according to Eq. (2.18) requires interpolations 
in three dimensions. We would like to reiterate here, that for the solution of the Faddeev 
equation, Eq. ( j2.14j ), only two-dimensional interpolations are required. To carry out the 
3D interpolation we use standard B-splines as given in ||22|| . In case of the wave function 
B-splines turned out to be slightly more accurate on the given grids compared to the Cubic 
Hermite splines. Since the wave function is calculated only once, the computational speed 
is not an issue here. 

The absolute value of the total wave function q, x) calculated from the MT-V po- 
tential for a fixed angle x — 1 is displayed in Fig. 4. This wave function has a smooth shape 
and is significant in size only for p and q values smaller than 1.5 fin -1 . It is interesting to 
see that the shape of the wave function is quite similar to that of the Faddeev amplitude of 
Fig. 2. The function values are of course different. 

In Fig. 5 the partial wave projected amplitudes ^i(p, q) are displayed as contour plots 
for / = 0, 2, and 4. The left column contains the ^i(p, q) obtained from the solution of 
Eq. ( p,19| ), whereas the right column contains ^i(p,q) calculated directly in partial waves 
from Eq. (2.24). Again, both plots are nearly everywhere indistinguishable, underlining the 
accuracy and consistency in the numerical realization of both schemes. For I = 4 the contour 
plot for ^a(p, q) calculated in partial waves show small irregularities for p values larger than 
5 fm -1 and q close to 0. These irregularities are missing in the projections obtained from the 
full 3D solution, but again as pointed out above the 3D solution is not correct there either. 

In Fig. 6 we show the wave function \l/(p, q, x) calculated from the MT-IV potential. A 
comparison with Fig. 4 shows that for the purely attractive potential the significant pieces 
of the wave function are more shifted to slightly larger momenta p and q. This is consistent 
with the larger binding energy given by the MT-IV potential, which causes the 3N system 
to be tighter bound when viewed in configuration space and thus giving a more extended 
structure in momentum space. Due to the lack of repulsion this wave function does not 
change sign. 



B. Momentum Distributions and Expectation Values 

In applications one does not access the wave function directly but rather only certain 
matrix elements thereof. One example, which can be 'measured' at least approximately in 
electron scattering is the momentum distribution, the probability to find a nucleon with 
momentum q in the nucleus. In our case this is given by 



10 



OO 1 

n(q) = 2nq 2 J dpp 2 J dx^ 2 (p, q, x) 

o -1 

1 °° 

= 4^ E <?' / <W ^ 2 (P, q, x). (4.4) 
? o 

It is of interest to see how the sum of the partial wave form approaches the expression 
evaluated directly from ^(p, q, x). In Fig. 7 we display the momentum distribution obtained 
from \l/(p, q, x) together with the partial wave sums corresponding to different /. The small 
q values are dominated by the / = part. However, for a correct representation of n(q) for 
larger q larger I values are required. While for q < 3 fm _1 I = and 2 are sufficient, at 
q ~ 10 fm _1 one needs I values as large as 8. 

Another interesting test for the numerical accuracy is the comparison of a complementary 
quantity, namely 

OO 1 

h{p) = 2n p 2 J dqq 2 J dx^ 2 (p,q,x) 
o -1 

1 °° 

= / dqq 2 ^ 2 (p,q,x). (4.5) 

1 o 

This quantity does not seem to be easily accessible experimentally. The momentum distri- 
bution n(p) as obtained from ty(p, q, x) together with its partial wave sums of different I is 
shown in Fig. 8. Here higher / values are mostly needed to fill in the dip at 2 fm _1 . 

A less detailed test for the quality of the 3N wave function is the evaluation of the 
expectation value (^/\H\^>) = (H) and compare this value to the previously calculated 
eigenvalues. Explicitly we evaluate the following expression 



= (^\H \y) + (V\V\V) = 3 (if)\H \^f) + 3 (4.6) 



where 



\H \V) = 8tv 2 J dp J dq 



OO OO _ „ 1 

2 „2 



1 ,2 i ^ 2 



— v + 

_m 4m 

-1 



dxp q ip(p, q, x) ^(p, q, x) (4.7) 



and 



CO CO OO 11 



(tflVtl*) = 8tt 2 J dpp 2 J dqq 2 J dp' p' 2 J dx J dx 

-1-1 

x ^(p, q, x) v\(p,p, x, x') ^(p', q, x). (4.8) 

Here vi(p,p',x,x') is the expression for the potential containing the integration over the 
azimuthal angle f (We use the same notation as given in Ref. 



(q r , q, x', x) = I dip V(q', q, x'x + y/l — x' 2 \/l — x 2 cos if). (4.9) 



o 
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In the case of the Malfliet-Tjon potential this integral can be performed analytically. 

The values of (H), (H ), (V) are given in Table V for both potentials MT-V and MT-IV 
calculated in both schemes. One can see that the energy expectation values and eigenvalues 
E t agree with high accuracy within each scheme as well as between the schemes. 



V. VERIFICATION OF THE 3N SCHRODINGER EQUATION 



In the previous sections we displayed properties of the 3N wave function, which we 
obtained from the Faddeev component by evaluating Eq. (2.18). We showed that the expec- 
tation value (H) evaluated with this wave function deviates from the calculated eigenvalue 
by less than 0.1%. However, since the expectation value is obtained by integrating over the 
wave function, this accuracy test gives information about the overall quality of the wave 
function. A more stringent test for the quality of ^(p, q, x) at each point in the p — q — x 
space is to determine the accuracy with which \l/(p, q, x) fulfills the 3N Schrodinger equation. 
To the best of our knowledge, such a test has never been carried out. 

The 3N Schrodinger equation for the bound state of 3 identical particles is given in 
Eq. (p.l|). Using the three different sets of Jacobi momenta in momentum space we obtain 



E t ^(p 1 ,q 1 ,x 1 ) 



-pi + 7— <?i 

m 4m 



#(pi,gi,zi) 

+ / dp'iPl J dx' 1 v 1 {p 1 ,p\,x 1 ,x' 1 )^{p' 1 ,qi,x' 1 ) 

-1 
00 1 



1 



h 2) 

-1 
00 1 

+ / dPsPs J dx 3 v 3 (p 3 ,p 3 ,x 3 ,x 3 )^(p 3 ,q 3 ,x 3 ). (5.1) 
-1 

Here Xi, 2^), i — 1,2,3 contains the integration over the azimuthal angle (p in the 

different subsystems % (Eq. ( |4.9| )). As an aside, in Ref. [|T^] this equation was used to obtain 
the eigenvalue E t and the 3N wave function. We insert ^>{p,q,x) obtained from Eq. (2.18) 
and verify the equivalence of both sides of the Eq. (5.1). For the numerical evaluation of the 
integrals in Eq. (5.1) the z direction for each integral is chosen separately, so that always q« 
points in the direction of z. Then the integration vectors p£ are chosen with respect to q« 
being the z-axis. The momenta and q^, % = 2, 3 have to be expressed as functions of pi 
and qi. From the definitions of the momentum vectors in Eq. (2.8) we obtain 

12 = pi - |qi 

P2 = -fPl - fqi 
Q3 = -Pi + |qi 

Pa = -|pi + fqi- (5.2) 
The corresponding angles are defined as Xj = pj • q^, where i = 2, 3. This leads to 
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%2 ~ 



y/Pi + ??i + 3xiP!q iy Jpl + \q\ - xmqi 

x 3 = pl-b!-xmqi (53) 

y^i + f<?i - 3xipigi^/p? + Jg? + rripigi 

Since the functional form of the integrals in Eq. (5.1) is the same, we only evaluate the 
first, which we can view as function of the momenta given in the coordinate system 1 

oo 1 

*(Pi,5i,Zi) =/ d P'iPi J dx[v 1 (p 1 ,p' 1 ,x 1 ,x' 1 )^(p' 1 ,q 1 ,x' 1 ). (5.4) 
o -1 

The two remaining integrals in Eq. (5.1) are then obtained by interpolating on x(p«> ^j)) 
i = 2,3 using the coordinate transformations given in Eqs. (5.2) and (5.3). 

In Fig. 9 we display the relative error A = \ (E^(p, q, x)—H^(p ) q, x))\/\E^(p ) q,x)\x 100 
for x = — 1 for the wave function obtained from the MT-IV potential. The graph shows 
that the 3N Schrodinger equation is fulfilled with a numerical accuracy better than 1% for 
momenta p and q up to 6 fm -1 . This is a typical range of momenta which enters calculations 
of matrix elements, e.g., for electron scattering. When p exceeds values of 10 fm" 1 the relative 
error becomes larger than 5% for almost all values of q. This behavior can be explained 
by inspecting the different Faddeev amplitudes separately as they are given in Eq. (2.19). 
If one considers the value x = — 1 as shown in Fig. 9, one sees that for p = 9.9 fm" 1 and 
q = 10.3 fm" 1 the Faddeev amplitude ip2 is needed at a value q2 = 15.05 fm" 1 , which is 
outside the g-range the Faddeev equation was solved in this particular case. Instead of 
extrapolating, we set ip2(p,q,x) = 0, which introduces a small discontinuity in the wave 
function. Though not visible in integrations over the wave function, this leads to a larger 
error when considering the pointwise accuracy of the 3N Schrodinger equation. Of course, 
this larger error in q, x) is completely irrelevant in all practical cases, since for those 
large momenta p and q the wave function ty(p,q,x) drops by about 8 orders of magnitude 
compared to its value at the origin. In the relevant momentum regions the 3N Schrodinger 
equation is fulfilled with very good accuracy by our three dimensional wave function. This 
statement holds for all x values. 

For the above study we chose the wave function obtained from the MT-IV potential, 
since this wave function does not have any node lines. It is quite obvious that our definition 
of the relative error would give a large error at the locations where ty(p,q,x) approaches 
zero. In fact, when calculating the relative error with the wave function obtained from the 
MT-V potential, we can identify the node lines of the wave function quite clearly. 



VI. SUMMARY 

An alternative approach to state-of-the-art three nucleon bound state calculations, which 
are based on solving the Faddeev equations in a partial wave truncated basis, is to work 
directly with momentum vector variables. We formulate the Faddeev equations for identical 
particles as function of vector Jacobi momenta, specifically the magnitudes of the momenta 
and the angle between them, and demonstrate their numerical feasibility and the accuracy 
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of their solutions. As two-body force we concentrated on a superposition of an attractive 
and repulsive Yukawa interaction, which is typical for nuclear physics, as well as on an 
attractive Yukawa interaction. The corresponding two-body t-matrix, which enters the 
Faddeev equations was also calculated as function of vector momenta. We neglected spin 
degrees of freedom in all our calculations. 

As first test for the numerical accuracy of the solution of the Faddeev equation as func- 
tion of vector variables, which is a three dimensional integral equation, we determined the 
energy eigenvalue of the bound system and compared our result with the one obtained in a 
traditional Faddeev calculation carried out in a partial wave truncated basis. We achieved 
excellent agreement between the two approaches as well as excellent agreement with cal- 
culations in the literature. We also found that the three dimensional Faddeev amplitude 
is nearly independent of the angle between the two Jacobi momenta. This is of course in 
agreement with the insights gained from the approach in the partial wave scheme. 

From the Faddeev amplitude we obtained the 3N wave function. We found here also 
that the dependence on the angle between the two Jacobi momenta is quite weak. We then 
performed a partial wave decomposition of our three dimensional solutions for the Faddeev 
amplitude and the 3N wave function and compared the so obtained partial wave amplitudes 
with the ones directly calculated. Again we found excellent agreement between the two 
approaches. In order to further probe the quality of our wave function we calculated the 
momentum distributions n(q) and n{p) and demonstrated that partial wave contributions 
up to I — 12 are necessary to build up the distribution n(q) for large momentum transfers 
q. In a similar vain we showed that higher partial waves are needed to fill in the sharp dip 
in n(p), which is obtained in a pure s-wave calculation. 

For a stringent test of the three dimensional wave function we inserted it into the 3N 
Schrodinger equation and evaluated the accuracy with which the eigenvalue equation is 
fulfilled throughout the entire space where the solution is defined. We found that within the 
physical relevant momentum region, namely p and q less than 10 fm _1 , the 3N Schrodinger 
equation is fulfilled with high accuracy by our numerical solution of the Faddeev equation. 

Summarizing we can state that the three dimensional Faddeev equation for a bound state 
can be handled in a straightforward and numerically reliable fashion. Once supplemented 
by spin degrees of freedom, this approach will most likely be more easily implemented than 
the traditional partial wave based method. State-of-the-art bound state calculations with 
realistic nuclear force models typically require at least 34 channels. The incorporation of 3N 
forces will most likely also be less cumbersome in a three dimensional approach. 
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APPENDIX A: EXPLICIT REPRESENTATION OF THE PARTIAL WAVE 
PROJECTED PERMUTATION OPERATOR 



For the convenience of the reader we provide the quantities G and G related to the 
permutation operators for general total angular momentum L. 

The function Gw(q,q',x') from Eq. (|2.15 ) is a combination of Legendre polynomials 

Pk(x): 



GVv(<7,g>)=E^) E E q'^q^g^r V2 - 

k Hl+H2=l Ul+U2=V 

For general L the geometrical coefficient is given as 



(Al) 



aSS^ - E k V/a/'A' ^ (2mi)!(2/X2)!(2 , i)!(2z/2)! 



x - 



IW\ 



1 A L g j \X'Lg'\ 1 V\ g' k j 

x C(/i 2 A g, 0) C(v 3 A' g\ 0) C(k v x g, 0) C(k ^ g', 0). 



(A2) 



Here I and A are the relative orbital angular momenta related to p and q. We also use the 
notation I = 21 + 1. 

In our context we have L = 0, which leads to I = A. Then g reduces to 



kil\V\p,2V2 

9iv 



IV 



l\l'\ 



\| fa h \| (2^ 1 )!(2/x 2 )!(2 I / 1 )!(2z/ 2 )! 

Y 1+Ul C(fi 2 I /xi, 0) C(u 2 I' vx, 0) C{fn v x k, 0) 2 (A3) 



V2+V2 



and Gix,i'\' to G w ._ 

The quantity G(p,q,x) occurs in Eq. (|2.24|) . In the case of general L the function G 
reads 

G IX Mp^)=Y, P ^) E E p^q^~g^r V2 ( A4 ) 

k 1-11+1-12=1' Ul+U2 = \' 
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with the geometrical factor 



9i\,vx 



= E k Vi'A'^^j (2/il)!(2/i2)!(2z/l)!(2z/2)! 

* „ V { 



g' g l 

l A k 



x C(//i z/i 0) C(// 2 ^2 0) C(</ fc Z, 0) C(g' fc A, 0). 
Again, for L = g reduces to 



9iv 



E*H 



Z'!A'! 



V 2 ^ \l A (2/i 1 )!(2/i 2 )!(2z/ 1 )!(2z/ 2 )! 



x 



(!) 



f /i! fX 2 I' 

\ v 2 v x g 
and correspondingly G^i'y to 



C(// 1 !/! 3, 0) C(/i 2 z/ 2 <?, 0) C(g k I, Of 
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TABLES 



TABLE I. Parameters of the Malfliet-Tjon type potentials. As conversion factor we use units 
such that he = 197.3 MeV fm = 1. We also use h 2 /m = 41.470 MeV fm 2 . 





V A [MeV fm] 


fiA [fm -1 ] 


Vr [MeV fm] 


MR [ fm *] 


MT-V 


-570.3316 


1.550 


1438.4812 


3.11 


MT-IV 


-65.1090 


0.633 







TABLE II. The calculated eigenvalue E t of the Faddeev equation as function of the number of 
grid points NP1, NP2, NQ1, NQ2, and NX chosen for the solution of Eq. ( gig ). The calculations 
are based on the potentials MT-V and MT-IV. 



NP1 


NP2 


NQ1 


NQ2 


NX 


Et [MeV] (MT-V) 


E t [MeV] (MT-IV) 


20 


12 


20 


12 


40 


-7.74387 


-25.0416 


24 


24 


24 


16 


40 


-7.74000 


-25.0453 


32 


24 


24 


16 


40 


-7.73761 


-25.0485 


32 


24 


32 


16 


40 


-7.73761 


-25.0483 


32 


24 


32 


32 


40 


-7.73761 


-25.0483 


32 


32 


32 


32 


40 


-7.73761 


-25.0483 


48 


32 


48 


32 


40 


-7.73666 


-25.0499 


64 


32 


48 


32 


40 


-7.73666 


-25.0502 


64 


32 


64 


32 


40 


-7.73650 


-25.0502 


72 


36 


64 


32 


40 


-7.73650 


-25.0499 



TABLE III. The calculated eigenvalue Et of the Faddeev equation as function of the number of 
partial waves employed. A fixed set of grid points is chosen corresponding to NP1 = 64, NP2 = 32, 
NQ1 = 64, NQ2 = 32, NX = 16. The partial wave i-matrix is calculated on a larger grid of 
160 x 160 grid points. The calculations are based on the potentials MT-V and MT-IV. 



1 


E t [MeV] (MT-V) 


E t [MeV] (MT-IV) 





-7.53975 


-24.8616 


2 


-7.71470 


-25.0465 


4 


-7.73383 


-25.0552 


6 


-7.73613 


-25.0562 


8 


-7.73649 


-25.0564 


10 


-7.73656 


-25.0565 


12 


-7.73658 


-25.0565 
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TABLE IV. The relative contributions of the partial wave Faddeev amplitudes Fi for each 
partial wave. We give also the corresponding values W\ for the wave function. The calculations 
are based on the potentials MT-V and MT-IV. 



7 

I 


Ft 


MT-V 


TT7" 

Wi 




F 


MT-IV 


TT7" 

Wi 







99.9951 






99.0851 




99.9940 






99.2121 




2 


0.4823 • 


10- 2 




0.7482 




0.5923 • 


io- 2 




0.7318 




4 


0.7988 • 


io- 4 




0.1159 




0.3128 • 


io- 4 




0.4682 • 


nr 1 


6 


0.2278 • 


10~ 5 




0.3305 • 


io- 1 


0.8813 • 


10~ 6 




0.6743 • 


10~ 2 


8 


0.1110 • 


IO" 6 




0.1088 • 


io- 1 


0.6056 • 


io- 7 




0.1520 • 


io- 2 


10 


0.8677 • 


io- 8 




0.3939 • 


io- 2 


0.7156 • 


io- 8 




0.4528 • 


IO" 3 


12 


0.1006 • 


io- 8 




0.1545 • 


io- 2 


0.1147 • 


io- 8 




0.1627- 


IO" 3 



TABLE V. The expectation value (H), (Hq), and (V) calculated for the potentials MT-V and 



MT-IV within the three dimensional scheme (3D) and the partial wave scheme (PW). 





{H ) [MeV] 


(V) [MeV] 


(H) [MeV] 


(E t - (H)) [MeV] 


MT-V 


3D 


29.77706 


-37.51340 


-7.73634 


0.00011 




PW 


29.7776 


-37.5139 


-7.73634 


0.00024 


MT-IV 


3D 


77.2055 


-102.2550 


-25.0495 


0.0003 




PW 


77.2653 


-102.3207 


-25.0554 


0.0011 
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FIGURES 



FIG. 1. The angular dependence for the symmetrized off-shell t-matrix, t s (p,p' ,x, E), is dis- 
played for p = 0.75 fm _1 and E = —50 MeV. 

FIG. 2. The magnitude of the Faddeev component tp(p,q,x) for x = 1 calculated from the 
MT-V potential. 

FIG. 3. The partial wave projected Faddeev components ipi(p,q) calculated from the MT-V 
potential are shown for I = 0, 2, 4. The left column of graphs represents the amplitudes obtained 
from the three dimensional calculations. Here (a) stands for the I = amplitude, (b) the one for 
/ = 2, and (c) the one for I = 4. The amplitudes in the right column are obtained as solutions of 
the partial wave Faddeev equation, (d) represent the I = amplitude, (e) the one for I = 2, and 
(f) the one for I = 4. 

FIG. 4. The magnitude of the 3N bound state wave function fy(p,q,x) for x = 1 calculated 
from the MT-V potential. 

FIG. 5. The partial wave projected bound state wave functions ^i(p,q) calculated from the 
MT-V potential are show for i = 0, 2, 4. The left column of graphs represents the amplitudes 
obtained from the three dimensional calculations. Here (a) stands for the I = amplitude, (b) 
the one for I = 2, and (c) the one for / = 4. The amplitudes in the right column are obtained as 
solutions of the partial wave Faddeev equation, (d) represent the / = amplitude, (e) the one for 
I = 2, and (f) the one for I = 4. 

FIG. 6. The magnitude of the 3N bound state wave function fy(p,q,x) for x = 1 calculated 
from the MT-IV potential. 

FIG. 7. The momentum distribution n(q) is shown in (a) for the MT-V potential as calculated 
in the two different schemes. The solid line refers shows the result from the three dimensional 
calculation whereas the dashed line gives the result from the partial wave sum up to I = 12. The 
monotonic increasing partial wave sums for / = 0, 2, 4, 6, 8, 10, and 12 are shown in (b). 

FIG. 8. The momentum distribution h(p) is shown in (a) for the MT-V potential as calculated 
in the two different schemes. The solid line refers shows the result from the three dimensional 
calculation whereas the dashed line gives the result from the partial wave sum up to I = 12. The 
monotonic increasing partial wave sums for / = 0, 2, 4, 6, 8, 10, and 12 are shown in (b). 

FIG. 9. The relative error A as defined in the text as function of the momenta p and q for 
the fixed angle x = — 1 represents a measure for the accuracy with which the 3N Schrodinger 
equation is fulfilled by the three dimensional wave function obtained from the MT-IV potential. 
The numbers at the contour lines give the relative error in %. 
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Fig. 1 
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Fig. 6 
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Fig. 9 
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